### Perform spatial cross-correlation test and visualization
### Stationary grid interpolation method (observations in different geo points)

# load libraries
library("spatialEco")
library("sp")
library("spdep")
library("gstat")
library("ggplot2")
library("tidyverse")
library("fields")
library("automap")
library("raster")

# read input point data
row1 <- read.csv("input/data-anchovy-lamp-2008-june.csv")
row2 <- read.csv("input/data-mnem-2008-june.csv")

# make coordinate-driven object
obj1 <- row1
obj2 <- row2
coordinates(obj1) <- ~x+y
coordinates(obj2) <- ~x+y

# perform spatial raster interpolation to get continious spatial function
# interpolation: by IWD, span = 0.5 (default recommended)

interp.box1 <- c(
  "xmin" = min(row1$x),
  "xmax" = max(row1$x),
  "ymin" = min(row1$y),
  "ymax" = max(row1$y)
)

interp.box2 <- c(
  "xmin" = min(row2$x),
  "xmax" = max(row2$x),
  "ymin" = min(row2$y),
  "ymax" = max(row2$y)
)

# calc step to find 100x100 box grid
interp.step1  <- c(
  "x" = as.numeric(interp.box1["xmax"] - interp.box1["xmin"])/100, ## change 50 to any other grid step, 50 = 50 * 50 cells
  "y" = as.numeric(interp.box1["ymax"] - interp.box1["ymin"])/100
)

# calc step to find 100x100 box grid
interp.step2  <- c(
  "x" = as.numeric(interp.box2["xmax"] - interp.box2["xmin"])/100, ## change 50 to any other grid step, 50 = 50 * 50 cells
  "y" = as.numeric(interp.box2["ymax"] - interp.box2["ymin"])/100
)

# reproduce grid 
sf.grid1 <- st_as_sf(row1, coords = c("x", "y"))
grid1 <- sf.grid1 %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_make_grid(
    cellsize = c(interp.step1["x"], interp.step1["y"]),
    what = "centers"
  ) %>%
  st_as_sf() %>%
  cbind(., st_coordinates(.)) %>%
  st_drop_geometry() %>%
  mutate(row1 = 0)

grid.raster1 <- grid1 %>%
  rasterFromXYZ()

sf.grid2 <- st_as_sf(row2, coords = c("x", "y"))
grid2 <- sf.grid2 %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_make_grid(
    cellsize = c(interp.step2["x"], interp.step2["y"]),
    what = "centers"
  ) %>%
  st_as_sf() %>%
  cbind(., st_coordinates(.)) %>%
  st_drop_geometry() %>%
  mutate(row2 = 0)

grid.raster2 <- grid1 %>%
  rasterFromXYZ()

# fit1 <- gstat(
#   formula = row1 ~ 1,
#   data = as(sf.grid1, "Spatial"),
#   nmax = 10, nmin = 3,
#   set = list(idp = 0.5)
# )
# idw1 <- interpolate(grid.raster1, fit1)
# plot(idw1, main = "IDW interpolation (w = 0.5, nmin=3, nmax=10)")
# 
# idw1.rms <- sqrt(mean(gstat.cv(fit1, debug.level = 0, random = F)$residual^2))
# 
# fit2 <- gstat(
#   formula = row2 ~ 1,
#   data = as(sf.grid2, "Spatial"),
#   nmax = 10, nmin = 3,
#   set = list(idp = 0.5)
# )
# idw2 <- interpolate(grid.raster2, fit2)
# plot(idw2, main = "IDW interpolation (w = 0.5, nmin=3, nmax=10)")
# idw2.rms <- sqrt(mean(gstat.cv(fit1, debug.level = 0, random = F)$residual^2))

## make stationary grid by minimum observed extent
xmin <- max(c(interp.box1["xmin"], interp.box2["xmin"]))
xmax <- min(c(interp.box1["xmax"], interp.box2["xmax"]))

ymin <- max(c(interp.box1["ymin"], interp.box2["ymin"]))
ymax <- min(c(interp.box1["ymax"], interp.box2["ymax"]))

# define point stationary grid steps
steps.x <- 40
steps.y <- 70

x.seq <- seq(xmin, xmax, by = (xmax - xmin)/steps.x)
y.seq <- seq(ymin, ymax, by = (ymax - ymin)/steps.y)

x.vec <- c()
y.vec <- c()

for (xi in x.seq) {
  for (yi in y.seq) {
    x.vec <- c(x.vec, xi)
    y.vec <- c(y.vec, yi)
  }
}

# make grid as data.frame for prediction
pts_grid.df <- data.frame(x = x.vec, y = y.vec)
pts_grid.sf <- pts_grid.df
coordinates(pts_grid.sf) <- ~x+y
#spplot(pts_grid.sf)

# interpolate by grid
fit1 <- autoKrige(
  formula = row1 ~ 1,
  input_data = as(sf.grid1, "Spatial"),
  new_data = pts_grid.sf
)

fit1.out <- fit1 %>%
  .$krige_output %>%
  as.data.frame() %>%
  dplyr::select(x = x, y = y, row1 = var1.pred, se = var1.stdev)

# iterp1 <- rasterFromXYZ(fit1)
# plot(iterp1, main = "Auto Kriging interpolation")

fit2 <- autoKrige(
  formula = row2 ~ 1,
  input_data = as(sf.grid2, "Spatial"),
  new_data = pts_grid.sf
)

fit2.out <- fit2 %>%
  .$krige_output %>%
  as.data.frame() %>%
  dplyr::select(x = x, y = y, row2 = var1.pred, se = var1.stdev)

iterp2 <- rasterFromXYZ(fit2)
plot(iterp2, main = "Auto Kriging interpolation")

## make prediction by interpolation layers by stationary grid
# pred1 <- predict(iterp2, pts_grid.sf)
# pred2 <- predict(fit2, pts_grid.sf)
# 
# data <- data.frame(x = pts_grid.df$x, y = pts_grid.df$y, row1 = pred1$var1.pred, row2 = pred2$var1.pred)
# geodata <- data
# coordinates(geodata) <- ~x+y

# data <- data.table::data.table(x = numeric(), y = numeric(), row1 = numeric(), row2 = numeric(), se1 = numeric(), se2 = numeric())
# data <- as.data.frame(data)

data <- data.frame(x = fit1.out$x, y = fit1.out$y, row1 = fit1.out$row1, row2 = fit2.out$row2, se1 = fit1.out$se, se2 = fit2.out$se)

# find "bad rows" interpolation if SE > median(se)
# and drop bad points from assessment
badrows <- c()
se1.q <- as.numeric(quantile(data$se1)[3])
se2.q <- as.numeric(quantile(data$se2)[3])
for (ri in 1:nrow(data)) {
  if (data[ri, "se1"] >= se1.q || data[ri, "se2"] > se2.q)
    badrows <- c(badrows, ri)
}

data.clear <- data[-badrows,]
geodata <- data.clear
coordinates(geodata) <- ~x+y

## perform cross-correlation test
w.euclid <- spDists(geodata)
fit.cor <- crossCorrelation(geodata$row1, geodata$row2, 
                            w = w.euclid, 
                            dist.function = "none",
                            k = 99,
                            scale.partial = T,
                            type = "LSCI")

#data$LSCI <- fit.cor$SCI[, "lsci.xy"]
geodata$LSCI <- fit.cor$SCI[, "lsci.xy"]

# vizualize
spplot(geodata, "LSCI", main = "Spatial CC (interpolated)")
morans.plot(x = geodata$row1, y = geodata$row2, coords = coordinates(geodata), scale.morans = T)
